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ABSTRACT 

'■>  This  paper  deals  with  the  numerical  solution  of  boundary  value  problems 
of  ordinary  differential  equations  posed  on  infinite  intervals.  The  solution 
of  these  problems  proceeds  in  two  steps.  The  first  is  to  cut  the  infinite 
interval  at  a  finite,  large  enough  point  and  to  insert  additional,  so  called 
asymptotic  boundary  conditions  at  the  far  (right)  end;  the  second  is  to  solve 
the  resulting  two  point  boundary  value  problem  by  a  numerical  method,  for 
example  a  difference  scheme.  In  this  paper  the  Box-scheme  is  investigated. 
Numerical  problems  arise,  because  standard  algorithms  use  too  many  grid  points 
as  the  length  of  the  interval  increases.  An  'asymptotic'  a  priori  mesh  size 
sequence  which  increases  exponentially,  and  which  therefore  only  employs  a 
reasonable  number  of  meshpoints,  is  developed.  ,  investigating  the 
conditioning  of  the  (linearized)  Box-scheme,  we  find  that  the  solutions  can  be 
obtained  safely  by  the  Newton  procedure  when  partial  pivoting  is  employed. 
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SIGNIFICANCE  AND  EXPLANATION 


This  paper  deals  with  the  numerical  solution  of  boundary  value  problems 
of  ordinary  differential  equations  posed  on  infinite  intervals.  These 
problems  have  the  following  form.  We  look  for  the  solution  of  a  system  of 
ordinary  differential  equations  which  is  defined  on  the  interval  [1,®]  , 
which  fulfills  a  certain  continuity  requirement  at  infinity  and  some  boundary 
conditions  at  t  *  1. 

Such  problems  occur  in  laminar  flow  theory,  in  fluid  dynamical  stability 
theory,  and  in  quantum  mechanics. 

For  the  numerical  solution  we  proceed  as  follows.  First  we  cut  the 
infinite  interval  at  a  finite,  large  number  t  *  T  and  impose  additional,  so 
called  asymptotic  boundary  conditions  at  t  =  T  in  order  to  obtain  a  'finite' 
two  point  boundary  value  problem.  Then  we  solve  this  problem  with  a  finite 
difference  method.  The  difficulty  arising  frequently  is  that  the  number  of 
mesh  points  has  to  increase  rapidly  as  T  +  »  in  order  to  preserve  a  certain 
accuracy  when  standard  algorithms  are  used. 


In  this  paper  we  derive  a  sequence  of  mesh- sizes  which  increases 
exponentially.  The  amount  of  computing  necessary  is  kept  reasonably  low. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


THE  NUMERICAL  SOLUTION  OF  BOUNDARY  VALUE  PROBLEMS  ON  ’LONG’  INTERVALS 
Peter  A.  Markowich*  and  Christian  A.  Ringhofer** 

1.  Introduction. 

In  this  pacer  the  numerical  solution  of  boundary  value  problems  on  infinite  intervals 
of  the  form 

(1.1)  y'  “  t“f(t,y),  1  <  t  <  »  ,  •  a  >  0 

(1.2)  b ( y ( 1 ) )  =  0 

(1.3)  y  S  C(  ( 1,*»] )  i<“>  y  8  C(  ( 1 ,  «>) )  and  lira  y(t)  -  y( «)  exists 

n+1  n  n  k  t*m 

is  considered.  Here  f:R  *  R  ,  b:R  ♦  R  where  generally  k  <  n  holds  because  (1.3) 

furnishes  another  set  of  boundary  conditions,  f  fulfills  certain  continuity  properties  at 

3f 

infinity  which  will  be  defined  later.  We  assume  that  the  Jacobian  —  (oo,y(«))  has  no 

3y 

eigenvalue  on  the  imaginary  axis. 

For  a  >  -1  Equation  (1.1)  has  a  singularity  of  the  second  kind  of  rank 
a  +  1  at  t  •  «>  .  But  we  disregard  the  practically  unimportant  case  -1  <  a  <  0  for  the 
following. 

Problems  of  this  kind  often  occur  in  fluid  dynamics  (boundary  layer  theory) ,  quantum- 
mechanics  and  electronics.  For  applications  see  Markowich  ( 1980a, b),  de  Hoog  and  Weiss 
(1980a),  McLeod  (1969)  and  Schneider  (1979). 

For  the  numerical  solution  we  proceed  as  follows.  First  the  infinite  interval  is 
substituted  by  a  finite  but  large  interval  and  n-k  additional,  so  called  asymptotic 
boundary  conditions  which  reflect  the  asymptotic  behaviour  of  the  solution  y  ,  are  imposed 
at  the  right  (far)  endpoint  T  .  We  obtain  two  point  boundary  value  problems  of  the  form 

(1.4)  x'  -  t“f(t,x),  1  <  t<  T 

(1.5)  b(x( 1 ) )  =  0 

(1.6)  S(T)x(T)  =  y(T) . 

t 
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The  condition  (1.6)  has  to  be  dhosen  such  that 


(1.7)  'Y  -  *°  88  T*“ 

holds  and  its  construction  is  described  in  de  Hoog  and  Weiss  (1980b),  Lentini  and  Keller 
(1980)  and  Markowich  (1980b). 


The  two  point  boundary  value  problem  (1.4),  (1.5),  (1.6)  now  has  to  be  solved  by  an 

appropriate  numerical  method,  for  example  by  the  Box-scheme  which  has  the  form 

x.  -  x. 

M.8)  ^  “  titV2f<ti+  V,  '  V2(Xi+1  +  Xi)>'  i  «  0(1)  (N-1 ) 


(1.9)  b(xQ)  -  0 

(1.10)  S(T)x  -  y(T) 

"  hi 

where  tQ  ”  1  <  t^  <  •••  <  <  tjj  “  T  ,  ”  fci  +  ^i  '  ^i+  l/^  “  t^  +  —  holds. 

It  is  clear  that  the  mesh-size  selection  is,  especially  for  these  problems,  very 

important  since  the  amount  of  labor  will  be  very  large  for  long  intervals  and  bad  (too 

small)  mesh-size  choices.  Even  well  working  adaptive  codes,  which  assume  a  relation  of  the 

form 

(1.11)  max  h  /min  h  <  const 

i  i 

and  whose  convergence  estimates  are  formulated  in  terms  of  max  h  ,  would  emmploy  too 

i  1 

many  meshpoints  in  order  to  admit  a  given  bound  for  the  global  error.  Moreover  the  codes 
which  employ  adaptive  mesh  refinement  (see  Lentini  and  Pereyra  (1977))  solve  first  with  a 
coarse  grid  in  order  to  do  local  error  estimation.  However  if  T  is  large  even  a  coarse 
grid  implies  a  lot  of  computational  labor  and  is  therefore  not  suitable. 

In  this  paper  we  use  the  asymptotic  form  of  the  solution  of  (1.1),  (1.2),  (1.3)  in 
order  to  construct  an  asymptotic  a  priori  mesh.  The  term  asymptotic  has  to  be  understood 
in  the  following  sense.  In  an  interval  (1,y)  ,  where  y  solely  depends  on  the 
•infinite'  problem  either  constant  meshsizes  or  more  sophisticated  algorithms  like 
equidistributlng  meshes  (see  Lentini  and  Pereyra  (1977))  have  to  be  used  and  in  ty.T)  the 
asymptotic  mesh  is  employed. 
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It  turns  out  that  mesh sizes  which  increase  exponentially  can  be  used  since  our 
assumptions  guarantee  that  y(t)  ■»  y(«)  exponentially.  For  the  Box-scheme  it  will  be  shown 
that  the  number  of  grid  points  which  is  necessary  in  order  to  achieve  a  total  accuracy 

0(e)  (total  accuracy  refers  to  the  difference  between  the  'infinite'  solution  y(t^)  and 

-  Vs 

the  discrete  approximation  x^)  equals  0( e  2  )  .  For  this  a  suitable  T  -  T(e)  has  been 
taken.  Moreover  it  will  be  shown  that  the  Newton  procedure  for  (1.8),  (1.9),  (1.10)  with 
these  exponentially  increasing  stepsizes  converges  quadratically  from  a  domain  of  starting 

values  which  does  not  shrink  to  4  as  e  ♦  0  .  He  also  show  that  the  condition  number  of 

-Vs 

the  (linear)  Box-scheme  is  an  0(e  2)  so  that  the  linear  system  can  be  safely  solved  by 

partial  pivoting.  Therefore,  in  the  nonlinear  case,  the  Newton  procedure  can  be  safely 
applied. 

Of  course,  no  fully  implicit  difference  scheme  (like  the  implicit  Euler  scheme)  should 
be  used  for  the  integration  of  (1.4),  (1.5),  (1.6)  since  in  general  the  fundamental  matrix 
of  the  linearized  problem  (1.4)  contains  exponentially  increasing  columns  which  are  scaled 
down  by  the  boundary  condition  (1.6).  But  nevertheless  this  would  cause  instabilities 
during  the  integration  when  using  large  meshsizes.  Similar  instabilities  can  occur  in  the 
case  a  >  0  when  using  the  trapezoidal  rule. 

Higher  order  stable  methods  can  be  constructed  by  polynomial  collocation  at  Gaussian 
points  and  will  be  analyzed  in  a  subsequent  paper. 

Another  way  to  solve  problems  of  the  kind  (1.1),  (1.2),  (1.3)  is  to  transform  the 

_  g 

'infinite'  problem  by  a  transformation  t  -  s  ,  (5  >  0  to  the  interval  (0,1)  and  to 
employ  difference  methods  at  this  constant  interval.  Methods  of  this  kind  have  been 
investigated  by  de  Hoog  and  Weiss  (1979).  However  this  way  of  proceeding  has  the 
disadvantage  that  a  singular  problem  (the  right  hand  side  of  the  equation  is  not  defined 
in  s  -  0)  has  to  be  solved  and  therefore  the  obtained  convergence  estimates  are  not  very 
strong.  Another  disadvantage  is  that  many  physical  problems  are  actually  posed  on  an 
Infinite  interval  (for  example  in  boundary  layer  theory)  such  that  a  'direct'  solution  is 


desirable . 


We  remark  that  there  la  a  cloae  connection  to  lingular  perturbation  probleme  aince  the 
traneformation  s  ■  j-  ,  e  “  takea  (1.4),  (1.5),  (1.6)  into 

(1.12)  e3^1*'  (a)  -  (a+£)°f(—  ,  «(a)),  0  <  a  <  1  ,  a  >  0 

e 

(1.13)  b( z(0 ) )  -  0 

(1.14)  e(“  +  1  )*( 1 )  -  0. 

(Note  that  lim  ,z )  -  f(®,z)). 

£♦0  E 

The  already  developed  mesh-aize  sequences  for  singular  perturbation  problems  cannot  be 
applied  aince  the  linearization  of  the  right  hand  side  of  (1.12)  does  not  generally  have  a 
series  expansion  in  powers  of  £  uniformly  in  0  «  a  <  1.  (see  Ringhofer  (1981))  since 
for  most  practical  problems 

(1.15)  £<t,y)  "  l  A1(y)t"i,  t  ♦  - 

y  1-0 

holds. 

Recently  Ascher  and  Weiss  (1981)  came  up  with  a  meshsize  sequence  for  linear  conatant 

coefficient  singular  perturbation  problems  (a"  0)  which  is  (as  a  formula)  equivalent  to 

ours,  however  they  apply  it  in  the  layer  of  thickness  0(e)  which  corresponds  to  the 
interval  t1,Y)  in  our  long  interval  case.  Outside  the  layer  they  use  a  uniform 
(independent  of  e)  mesh  fine  enou  jh  to  approximate  the  solution  of  the  reduced  problem 
(e  -  0)  well.  The  difference  comes  from  the  fact  that  the  solution  of  the  singular 
perturbation  problem  decays  exponentially  to  the  solution  of  the  reduced  problem  within  the 
layer  while  in  the  infinite  Interval  case  exponential  convergence  holds,  for  a  *  1. 

Krelsa  (1975)  used  a  similar  approach  to  construct  meshes  for  singular  perturbation 
problems. 

This  paper  is  organized  as  follows.  Chapter  2  gives  a  short  summary  of  the  theory  of 
boundary  value  problems  on  infinite  intervals  and  their  'finite'  approximation,  in  Chapter 
3  stepsize  sequences  are  developed  for  scalar  initial  value  problems,  Chapters  4,5  deal 
with  linear  boundary  value  problems  and  in  Chapter  6  nonlinear  problems  are  dealt  with.  In 
Chapter  7  the  results  are  gathered  and  put  into  algorithmic  form. 
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2.  Boundary  Value  Problems  on  Infinite  Intervals  and  their  Approximation  by  ‘finite* 
Interval  Problems  -  A  Summary 

He  consider  boundary  value  problems  on  infinite  intervals  of  the  following  form 

(2.1)  y'  -  t°Mt)y  +  t“f(t),  1  <  t  <  »  ,  a  >  0 

(2.2)  By(1)  -  S 

(2.3)  y  e  C([1,«l) 

where  the  nxn-matrix  A  e  C(f1,<»l)  and  f  e  C([1,«])*  B  is  a  matrix  whose  rank  is  (in 
general)  less  than  n  since  (2.3)  furnishes  another  set  of  boundary  conditions. 

Let  us  first  consider  the  case  where  A(t)  =  A  .  A  shall  have  the  Jordan  form  J 
obtained  by 

(2.4)  A  =  E  J  E*1  . 

We  assume  that  J  has  the  block  form 


(2.5) 


J 

*-» 

r 


K 

}r 


where  the  r+  x  r+matrix  J+  has  only  eigenvalues  with  positive  real  parts  and  the 
r_  x  r  -  matrix  J~  has  only  eigenvalues  with  negative  real  parts.  Imaginary  eigenvalues 
will  be  excluded  for  the  following.  The  diagonal  projection  D+,  D_  are  defined  by 


(2.6) 


0 


(2.7) 


The  general  solution  of  (2.1)  (with  A(t)  =  A)  and  (2.3)  can  now  be  written  as 


(2.8) 


y(t) 


E$(t) 


?  +  E( HE* 1 f ) ( t ) . 


r 

(6C' 
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--“v-  -  •  •’ 


where 


(2.9)  4>{t)  =  expiry  ta+1) 

is  Che  fundamental  matrix  of  the  transformed  problem 

(2.10)  u'  -  taJu  +  t^E^ftt) 

where  Eu  ”  y  holds  and  (HE_1f)(t)  is  a  suitable  particular  solution  of  (2.10)  which  can 
be  taken  as 

(2.11)  ( Hg )  ( t )  -  <J>( t )  D+$  1(s)s0tg(s)ds  +  $(t)  /tD_4>  1  ( s)sag(s)ds 

“  Y 

for  some  y  >  1.  This  operator  has  been  analyzed  by  de  Hoog  and  Weiss  ( 1980a, b)  and 

Markowich  (1980a). 

H  has  the  following  properties 

(  2. 12)  (a)  HsCUy,-])  ♦  C([Y,»]> 

(2.121(b)  IHI.  ,  <  c 

LY<°“J 

where  1*1,  denotes  the  max-norm  on  (Y/“l  reap,  the  associated  operator-norm.  The 

l  Y»“l 

constant  c  is  independent  of  Y  • 


Markowich  (1980a)  has  shown  that  if 


(2.13) 


f(t)  -  F(t)exp(-  -fr  t"  ),  F  e  L  ([1,«]  )  n  C(  [1,->> 

CfT  1  00 


holds  with  U  >  X  ,  >0  where  X  ,  is  the  smallest  modulus  of  the  real  parts  of  the 

min  min 

eigenvalues  of  A  which  are  in  the  left  half  plane,  then 


(t)l  <  const. (»FI  +  ISO  »4>(t )  T°  I 

[I,*]  1^ 


,  ..  ,,<?■  .  ,,n.  r-'  “*')  „  i  min  .  orM. 

<  const. (IFI,.  ,  +  l?l)t  exp(-  ■  t  ) 

[I,*]  cr-l 

holds,  r  is  the  largest  dimension  of  a  Jordan  block  associated  with  an  eigenvalue  of  A 
with  real  part  "*min  • 

The  boundary  value  problem  (2.1),  (2.2),  (2.3)  with  A(t)  =  A  is  uniquely  soluble  for 


all  $  e  R  ”,  f  e  C( [1,»] )  if  and  only  if  the  r  x  r  matrix 
(2.15)  BE+(1)[  °  1  is  nonsingular. 


Here  B  is  assumed  to  be  a  r  x  n-matrix.  So  the  continuity  requirement  (2.3) 


furnishes  r+  linearly  independent  boundary  conditions. 
(2.15)  implies  that  1£l  <  const. 161  holds. 


The  variable  coefficient  case  A(t)  |  A  is  treated  by  a  perturbation  approach  (see  de 


Hoog  and  Weiss  (1980a, b)  and  Markowich  (1980a, b)).  A(«)  now  plays  the  role  of  A  .  We 
assume  that  A(  <•)  has  the  Jordan  form  J  given  by  (2.5).  Then  we  can  show  (see  de  Hoog 
and  Weiss  ( 1980a, b))  that 


(2.16) 

y(t)  » 

E^_(t)c  +  E*(E  f)(t),  sec 

where  (t) 

is  an 

n  x  r_-matrix  defined  by 

(2.17) 

<!>_(•) 

“  (I  -  H(E_1A(.)E  -  J))"V*)[  J 

]  e  C<  [Y,»l  ) 

r 


for  y  sufficiently  large.  For  t  e  [1,y]  i(i  tt)  can  be  continuously  extended.  E^i(E  ^ f ) 

is  a  suitable  particular  solution  of  (2.1)  .  The  boundary  value  problem  (2.1),  (2.2), 

r_ 

(2.3)  is  uniquely  soluble  for  every  6  8  R  ,  f  ec([1,*)|  if  and  only  if  the 
r  x  r  -matrix 

(2.18)  BEi|r  (1)  is  nonsingular. 

Markowich  (1980a)  has  proven  the  estimate 

(X  -6) 

(2.19)  ly(t)  I  <  const  (IFIt1  +  lgl)exp( - ^2 -  t0^  ) 

for  t  >  t  .  5  >  0  can  be  chosen  sufficiently  small  if  t  is  large. 

Now  we  briefly  consider  nonlinear  problems  of  the  form  (1.1),  (1.2),  (1.3). 

From  (1.1),  (1.3)  we  conclude  that 

(2.20)  f ( «,y( “) )  -  0 

has  to  hold.  We  assume  that  the  roots  y(w)  of  (2.20)  are  isolated  and  take  one 

*  * 

particular  root  y  («>)  for  the  following.  Moreover  f(t,y  («•))  shall  fulfill  (2.13). 
Defining 

C^(t,a)  =*  { (t ,y)  e  Rn+'|t  >  t  ,  ly  -  y*(»)l  <  k) 

* 

we  assume  that  f,  fy  e  Clip(CM1,y  (■)))  for  a  sufficiently  large  k  .  We  also  assume 
that  the  boundary  value  problems  (1.1),  (1,2),.  (1.3)  has  an  isolated  solution,  i.e.  the 


linearized  problem  is  nonsingular. 


Now  let  J  be  the  Jordan  form  of  f^(“,y  («•))  obtained  by 
f  (»,/(«>)  =  EJE_1 

and  let  J  fulfill  (2.5)  such  that  D,,  0  ,  are  defined  as  in  (2.6),  (2.7).  Then  X  . 

"■  mm 


* 

is  defined  for  J  as  above  and  $(t),  <J)_ ( t )  are  as  in  (2.9),  (2.17)  with  f^(t,y  (“)) 
substituted  for  A(t).  Markowich  (1980a)  showed  that 

(2.21)  ly(t)  -  y  (»)1  <  const.  l\|i_(t)  I  <  const.  exp( - ~^+l —  tC<+ 

The  isolatedness  of  y  now  implies  that 

(2.22)  (y(  1 ) ) i|>  (1) 

dy 

is  nonsingular.  More  information  on  the  analysis  of  these  problems  can  be  found  in  the 
above  cited  references  and  in  Lentini  and  Keller  (1980). 

We  want  to  approximate  the  1  infinite'  problems  (2.1),  (2.2),  (2.3)  by  'finite'  two 
point  boundary  value  problems  of  the  form 

(2.23)  x'  =  taA(t)x  +  taf(t),  1  <  t  <  T  ,  T  >>  1 

(2.24)  Bx ( 1 )  =  B 

(2.25)  S(T)x( t)  =  y(T)  . 

Since  (2.24)  is  a  boundary  condition  of  rank  r_  we  assume  the  S(T)  is  a 
r+  x  n-matrix.  The  question  that  arises  immediately  is  how  to  construct  an  asymptotic 
boundary  condition  S(T)  such  that 

(2.26)  ly  -  xl[1iT]  *0  as  T  ♦  - 

and  the  order  of  convergence  should  be  as  fast  as  possible. 

A  complete  theory  of  this  kind  can  be  found  in  de  Hoog  and  Weiss  (1980a)  and  Markowich 
(1980b)  and  therefore  we  only  give  excerpts  which  will  be  needed  in  the  sequel.  The  basic 
idea  is  that  the  boundary  condition  (2.25)  has  to  scale  down  all  solution  components  of 
(2.23),  which  do  not  decay  exponentially. 

We  assume  that  (2.18)  holds.  A  possible  choice  is 

(2.27)  S  =  SIT)  -  (I  ,  0]E-1,  y(T)  =  0  . 

+ 

It  ha3  been  shown  in  the  above  cited  references  that  this  boundary  condition  implies 
convergence  in  the  sense  of  (2.26)  and  that  for  general  y(T) 

(2.28)  ly  -  xl  <  const ISy(T )  -  y(T)  I 

holds.  In  general  the  admissibility  conditions  for  a  boundary  condition  S(T)  are 

(2.29)  1S(T )  I  <  const.,  T  *  » 


(2.30) 


)  1  I  <  const. , 


T  ♦ 


I  (s(T  )E 


then  (2.28)  holds  for  the  unique  solution  x  of  the  'finite'  problem  if  T  is 
sufficiently  large.  y(T)  =  0  is  a  natural  choice  for  linear  problems. 

If  f(t)  fulfills  (2.13)  an  estimate  for  the  order  of  convergence  is  given  by  the 
right  hand  side  of  (2.19).  Moreover  it  has  been  shown  that  the  choice  (2.27)  is  optimal  in 
the  sense  that  the  actual  order  of  convergence  exceeds  (2.19)  for  homogenous  problems. 

The  condition  (2.25)  with  S  fulfilling  (2.29),  (2.30)  and  y(T)  -  S(T)y  (»)  can  also 
be  used  for  nonlinear  problems  of  the  form  (1.1),  (1.2),  (1.3)  if  the  above  stated 
assumptions  on  f(t,y)  and  the  solution  y  hold.  (2.28)  still  holds  for  nonlinear 
problems . 
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3.  The  Scalar  Case 


In  this  chapter  we  treat  the  simplest  case,  namely  scalar  initial  value  problems.  ' 
aim  is  to  construct  step-size  sequences  for  the  Box-scheme  such  that  the  global  error  is 
less  than  a  prescribed  accuracy  regardless  of  the  length  of  the  interval  of  integration. 
These  step-size  sequences  will  be  used  for  the  general  boundary  value  problem  case. 

He  consider 


(3.1)  y’  -  -Xt°y  +  t“f(t)  ,  1  <  t  <  -  ,  a  >  0 


(3.2) 


yd)  »  y 


where  X  •  X^  +  iX^  may  vary  in  a  compact  subset  ft  of  {z  e  C|Rez  >0} 


(3.3) 


The  Box-  or  centered  Euler  scheme  for  (3.1),  (3.2)  has  the  form 

a 


yi+^i  X 


'  2  fci+  Wyi+1  +  yi)  +  *i+  V,  fi+  V,  '  1  >  0  »  y0  =  y 


where  for  h^  >  0 


'2  '2 
h. 


0  "  “i+1  “i  '  "i  '  ~i+ V2  i  2 


(3.4)  (a)  t„  -  1,  ttJ_,  -  t,  +  h,  ,  t 

holds. 

We  define 


t.  +  —  ,  i  >  Of  (b)  f 


i+  V2  "  f(ti+  V2 


(3.5)(a) 


m  1  -  |  h  t®  i. 

y  ( X,h)  =  n  - .  -  ,  n  <  m  ,  Y  -  1  ,  n  >  -1 

n'm  . .  X  a  n+1,n 

*  "  1  +  -  httt+  y2 


i-1 


and  for  a  sequence  of  complex  numbers  z  *  (as.  4/). 

J+  '2  J*1 


i-1  h  t  1/  z  1/ 

( 3.5 )  (b)  (H  ( X,t  ,h)z)  =  l  1  &  —2  Y  ,(X,h),  i  >  1+1 


i-I  1  +  2  h4tt+  y2 


.i-1 


and  (H_(X,tI,h)z)I  =  0  where  h  -  (h  )^_0  is  the  sequence  of  step-sizes. 


Using  these  definitions  the  solution  of  (3.4)  can  be  written  as 


(3.6) 


yi  “  Y0,i-1U'h)Y  +  (H-(X'Vh)f,i 


where  f  ■  ( f 


j+  v2  *j 


i-1 


has  been  set. 


The  local  discretization  error  t®+  y^  t1+  y^  of  the  difference  scheme  (3.3)  Is  defined 


(3.7) 


y<t.  )-y(t  )  . 

fci+  V,  Vj“  hT  +  2  tl+  V2y<ti+1)  +  y(ti)>  "  t?+  V.  fi  +  I/-'  1  ?  0 


72  if  72  n^  2  1+  V2 

where  y  is  the  solution  of  (3.1),  (3.2) 

The  global  error 

(3.8)  e^  «  yft^)  -  yA 
then  satisfies  the  difference  equation 

(3.9) 


i+  V2  i+  V 2 


1  ‘  2  t“+  V2  (°i-M  +  V  + 


•1  *  ■■  '2 
and  therefore  has  the  solution 


i+  V2  *i+  V2 


i  >  0  ,  eQ  «  0 


(3.10) 
with  f 


(H_U,t  ,h)t) 


.1-1 


( 1 j+  V2  i-o  • 

In  order  to  estimate  the  right  hand  side  of  (3.10)  we  need  the  following 


Lemma  3.1. 


Let  x^  for  x  -  n,n+1,.«.,m  be  complex  numbers  with  Rex^  >  0  and 


Imx 

_ _ K 

Rex 


<  x  »  const.  Then  setting 


i 

n 

i-1+1 


0  for  i  >  1 


I* 


(3.11) 


K*n  |l+x  I2  i»x+1  ,1+Xil 


1"X1I  ,  / - 2 

- <  V2  /l+x 


(3.12) 

holds. 

Proof: 


,Xx' 


m 

1  , 1  r 
■c*n  1 1  +x  I 

< 


K-1  I  1 —X  | 

n 

t=n 


'1+Xil 


<  V2 


1+x 


An  easy  calculation  gives 


Ixx' 

|Uxx': 


/ - 5  |1“*J 

<V2/1+x  (l.-pj^j- 


)  . 


Substitution  into  the  right  hand  sides  of  (3.11),  (3.12)  yields  telescoping  sums. 
Application  of  Lemma  3.1  immediately  yields 

Lemma  3.2.  Let  f  =  (f.  1.)^  2  •  Then  for  every  sequence  h  -  (h.)*  1  with  h.  >  0 
-  1+  /2  £=I  i  £»I  t 


(3.13) 


I  (H  (X/t  »h)f)  I  4  const.  max  (If..  ,,  Id  +  h«t?+  1/>1 

-  I  1  t-I(1)(i-1)  t+/2  2  1 


holds  for  i  >  I  uniformly  for  X  e  fl  • 
Proof : 

|(H_(X,tI,h)f)il  < 


4  const  max  ( I  € 0  1/ M  1  +  h  . t a  ^.) )  V  - —  —  •2~ 

„  1+  Vo  2  l  1+ y2  ...i.  v.  ta 


w(1)(i-1l  "  '2 

and  application  of  (3.11)  yields  (3.13). 


*-i  |1+1/2Vi+v,' 


2  Vl,Ha'h) 


|e  |  4  const.  max  ( •  ^.4. 1/1 (1  +  h»t?+ V  )]* 

1  Wdlli-ll  1  '2  11/2 


We  get  from  (3.10) 

(3.14) 

l*0(1)(i-1) 

For  the  following  we  assume  that 

(3.15)  f(t)  *  F(t)exp(-  t®*1)  ,  u  >  X, 
holds  where  F,F',F”  e  C([1,-))  n  L  ((1,»J). 


(3.16) 


R  straightforward  calculation  gives 

2r  1_ 

1  V2 


1 1,+  1/  I  <  const.  h2[— ly"*l  t  ,  +  Ully’l,  t  .]• 

U  x2  4  t“+V2  ‘Wl1  ‘Wl1 


Markowich  (1980a)  shows  that  (3.14)  implies 

X1  ™+1 

I  y(  t )  |  c  c^XK  IFl(1  ^  +  ly|)exp(-  —  t  ) 

where  c  (  X)  is  bounded  for  X  e  ft  •  Differentiating  (3.1)  and  using  (3.16)  yields 
'  2  X 

(3.17)  Ul+y2}  <C2(X,(  l  ,F<K),(1,-}  +  IvD^-^tf1)  •  tl>t(1)U) 

where  the  function  t^aexp(-  — t®*1)  takes  its  maximum  over  [1,«°]  at  t,,.(X). 

ortl  ( 1 ) 

c2(X)  and  tj^(X)  are  bounded  for  X  S  fl  .  We  get 

le  |  4  c(  l  IF(  K>  I  +  |y  |  Xtftl  max  h2  + 

1  k-0  ,1'*1  (1)  1*0  (1)1  4 


(3.18) 


max  [h2t2aexp(-  tf')(1  +  h  ,  )]  . 

l-(I+1)(1)(i-1)  *  Ctrl  X.  '2 


Here  t  <  t  4  t  holds,  where  t  *  max  t  (X)  and  c  is  independent  of 

•  1  I  f  1  '  I  l  I  /  Q  1  I  / 


x  e  a 


(3.18)  gives 


ej  <  2c (  l  IF<K,I 

K-0 


<*>■  +  I  vt )  (t3*1 

[1,.)  +  IY,Ut(1) 


max  h,  + 
t 

t-omi 


(3.19) 


max 

i-(I+1)(1)(i-1) 


[ttaexp('  7T  tfS-.<hJ.hK.hf3„ 


How  we  require  that  for  some  0  <  e  <  E 


(3.20) 


|e  |  <  4c[  l  IF(K,|.  .  +  )yi )e  ,  i  >  0 

<-0  1,1 
holda.  This  is  fulfilled  if  we  choose 


(3.2 1 ) (a) 


h,  <  , 

t-0(1)I  1  /t^| 


h(  X1#  Eit^),  t^  (  ^(1) 


and  for  t  >  I 


r ^p'IF^TT  tT1  fc(i)  <  fct  <  (ff  _  t(2>(E) 


( 3.21 )  (b)  h^  <  MA^e.t^)  -  / 


3/e  t(2)(E)  4  tt  «  t(3)(e) 


_  _2a 

o+3  /"  .  cr+3  , 
✓e  t  exp( 


X1  .ati. 

(at3)(at1)  ti  ’  '  ti>t(3)(E) 


where  t(3)(e)  >  t(2)(e)  is  the  root  of  -  t(3)*+1  >exPt3(^fiT  )  1 '  (Note  that 

/e 

t(2)(c)  -  t(3)(E)  for  a  -  0.) 

This  defines  upper  bounds  for  the  stepsizes  at  a  given  point  t^  depending  on  e  , 
the  bound  for  the  global  error,  and  on  A  *  Re  A  •  These  bounds  are  independent  of  the 
length  of  the  integration  interval  and  increase  (in  t  ”  t^)  exponentially. 

We  now  compute  the  number  of  steps  N(T, s)  which  is  necessary  for  integration  on  the 


interval  [1,TJ  if  h  »  h(A  ,e,t.)  is  chosen.  Therefore  we  write 
X  IX 


(3.22) 


where  I . 


N(T,e) 


l  1  +  l  1  +  I  1  +  I 


iSI, 


iSI, 


(n) 

set.  Obviously 


(0)  (1) 

U  e  ".It.  0(t,_.,t.  _ J >  where  t 


iei 


(2) 


iei 


(3) 


0’  1  (n ) '  (n+1) 


(0)  “  fc0 


1  and  t^4j  =  T  has  been 


(3.23) 


iei 


l  1  < 


(t 


(1) 


-iwG? 


(D 


(0) 


/e 
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( 3 . "4 ) 


hi  hi+1  .  ,hi+1.  ft(2)  dt  .  ...  1 

-r —  <  max  (-r — )  J  -  <  c(X)  - 

iei(1)  iei(1)  “i+1  1  i01{i>  1  fc(i)  s<*'e't> 


J  1  -  l  h 


i+ 1  2 

where  c(X)  is  bounded  in  fl  since  max  — —  <  e  holds.  The  same  estimate  is 

iei  hi 

easily  obtained  for  I  1  • 


iei 


(2) 


For  t^(e)  <  t^  <  T  we  get 


(3.25) 


-  X 

hi+1  *  2  *Xp((at3)(at1)  (2  _1)t(3))hi  ’ 


Since  tft?  >  -  in  —  holds  we  derive 

(3)  X,  e 


2a  (2af,-1) 


2a  1 


(3.26) 


,  .  ,  at  3  at  3  at  3  3 

h±+1  >2  e  hl  >  2  e  ht 


and  get  a  bound  K  «  £  1  from 

iSI(3) 


(3.27) 


4,+k-l 


T  > 


i-j 


I  \ 


where  I  »  +  K-t)  has  been  set.  Further  we  use 

1  1 

,  ,  ,  .  /at1\ott1  l.atl 

hj,  >t(3)(£l  >  )  (tne’ 


such  that 


(3.28) 


2g  _1  _1_  1 

a+3  3ii  <- ortl  ^ art  1  /■  1->a+1 


\ti  >  t2  E  1  "t 
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From  (3.27),  (3.28)  we  conclude 


The  solution  of  (3.34)  is  given  by 


(3.36) 


z^  “  Y^  +  (H+(m,t^,h)f  )^ 


where 


N-1  h  y  f  y 

(3.37)  (H  (u,tM,h)f)  -  -  l  —  *  /?-—  /?  Y.  ,  (dl,h) ,  i  <  N-1 


‘-1  1  +  f  Vi*  v2 


and  (H+(  u,tN,h)f )  -  0 


Since  the  increasing  components  are  scaled  down  by  the  asymptotic  boundary  conditions 
at  t  -  T  we  disregard  the  convergence  of  the  z^  to  z(t^)  but  we  prove  a  stability 
estimate  analogously  to  Lemma  3.2. 


Lemma  3.3.  Let  f  ■  (f,.  i .)'  .  Then  for  every  sequence  h  »  (h_)_  . 

1+  79  l"l  «  l  Zm  i 


(3.38)  I  (H+(oi,tN,h)f)1|  <  const.  max  1 1 f  *+  1/1 ( 1  +  ^  htfc 1+  V  ’ 1 


t-l(1)(N-1) 


holds  for  i  <  N  uniformly  for  u  6  S)  . 
Proof. 


I  (H+((o.tN,h)f)il  < 


lVv2l(1+Th4v2,,I  ,2  '*i,t-i<“'h>'' 

Application  of  Lemma  3.1,  (3.12),  yields: 


<  max  ( | 

Jt-i(  1 )  (N-1 )  *T  '2 


2  ‘-1 


|(H  (u,t  ,h)f)  |  <c(u)  max  If  1/  t  ( 1  +  -Mp  h,t“  ,  ) . 

N  1  i-i(IHN-l)  1/2  2  1  1  '2 


Finally  we  prove 


Lemma  3.4.  Assume  that  <  t^  <  t^fc)  and  that  h^  <  hfX^e.t^)  for  ^  >  0 


Then 

and  c  -  c(o>)  is  bounded  on  !1  . 


Y.  .  , ( u,h) |  <  exp(-c(t.  -  t  )) 
J  i 


Than 


Proof. 


holds. 

tit  m  til. 


Sines 


Let  z  -  r1  +  iz2  ,  z1  >  0  . 


|!=£|2 

1+z' 


4i. 


<  1  - 


h+*r 


<  axp(-4 


I 1+* I ‘ 


This  estimate  has  been  used  In  de  Hoog  and  tfelss  (1979].  Therefore  for 
I  +  lu2 


lYi,J-1<“'h> I 


<  exp(-2<i>1 


j-1 

I 


1-1 


h„t 


I  t*  'A 


11  +  a  «  c(“>  for  *i  4  fcj 


<  t2u> 


holds  we  get 


2  oi  j-1  j-1 

|Yl,j-iu'h"  <  «p<-  I  4  exP<-c  l  V  <  exP  (-c<tj  * 

XtmX 


4.  The  Case  hit)  S  A 


We  consider 

(4.1)  x'  -  t°Ax  +  t°f(t),  1  <  t  <  T  ,  a  >  0 

(4.2)  Bx( 1 )  -  B 

(4.3)  S(T)x(T)  -  -y(T> 

r_ 

where  A  fulfills  (2.4)/  (2.S),  B  Is  an  r  »  n-matrix,  g  e  R  and  f  fulfills 

(2.13)  with  p  >  X  .  ,  the  r^  x  n-matrix  S(T)  fulfills  (2.29)/  (2.30).  This  simple 

min  + 

case  shall  be  considered  as  a  model  for  problems  where  A  depends  on  t  . 

The  Box-scheme  has  the  form 


(4.4) 

(4.5) 

(4.6) 


X  —  X 

i+1  i  A  « 


2  *1+  v,<xiti +  v +  *:♦  v,fi+  vi  1  -  on)<»-1 > 


h^  x  if  -/2  ifi  l  if  y2  if  y2 

Bx0  -  B 


,N-1 


S(T)Xj,  ■  y(T) 

where  the  partition  A  «  ^o  ,fc1 '  *  * ’'Sl-I '  fcN^  fulfill*  (3.35)  and  h  -  ,  hA  >  0. 

A  fulfills  (2.4).  Now  we  employ  the  transformation 
(4.7)  xA  -  Buj^ 

and  get 


(4.8) 

(4.9) 

(4.10) 


“ifi  “  “i  J  a 


2  *:+  v2  (wv +  v2E_1fi+  v2 '  1  ■  o<1,(N-,, 

BEUp  "  B 


SfTJEUj,  -  Y(T)  . 

We  want  to  derive  an  existence  and  stability  theorem  for  (4.8),  (4.9),  (4.10).  As  de  Hoog 
and  Weiss  (1980a)  did  for  the  continuous  case  we  split  u^  into 


(4.11) 


4-  | 

ui  1 

>*♦ 

(E_1fitV2,+ 

ui  “ 

.  ui  J 

$ 

}r_ 

(b)  E"  fit  y2" 

(E'1fi+v2r 

and  get  employing  (3.6)  for  u^  resp  (3.36)  for  u^ 


<4.12) 


r  r+  ♦  i 

-  - 

Yi,N-1<J  <h> 

e  4, 

0 

r-  -  j- 

0 

Vi-i'-1  'h)a 

where  for  any  k  x  k  aatrix  P  whose  eigenvalues  have  positive  real,  part 
<4.13)l£  _<P, 

l-n 

ho Ida  and  tha  oparator  H  ia  dafinad  aa 


n,m<P'h>  ‘  “  ‘V  f  !VuV2  >'  “  »  n'  Vu,P'hl  ■  *'  n 


(4.141(a) 


0'  N' 


+  t 

(H+(J  ,tN,h)*  )t 

1 

+ 

z 

m 

9  Z  m  l 

^  (H_(-j",t0,h)r')i 

z 

and 


N- 1  r 


(4.141(b)  (H+(J+,tH,h)a+)1  -  -J1htYl!t-1,J+'h,(I+  I  httl  V2  r1t“+  V2  V  V; 

<4.141(0  (H.(-J,tI,h)a-)i  -  -  f  V*  V2*U  V2 


whara  *+  -  <  *t+  V^l-0  '  *’  *  (*i+ V,  >1-0  an<1  et+ V,  ®  C  +  '  *1+ V  «  C  "  hM  b88n  88t- 


Haro  (H+(J+,t1|,h)z+)s  »  0  and  <H_(-j”,tI,h)l”)I  ■  0  hold.  These  definitions  make  sense 
because  (I  +  tJ+)  (I  -  tJ  )  1  exist  for  t  >  0  . 

In  order  to  get  bounds  for  the  defined  difference  operators  we  use  the  following  well 
known  representation  of  a  matrix  function 


(4.15) 


*<P>  -  /  V*  X)  <  AI-P)"1d* 


where  tha  contour  encloses  all  eigenvalues  of  P  .  f  is  assumed  to  be  analytic. 

If  all  eigenvalues  of  P  have  positive  real  parts  we  get 

(4.16)  Y*  (P.h)  -  -~r  !  T  mU,h)<U.-P)'1dJ 

n»w  zitx  p  n,m  ^ 

where  Y_  is  defined  in  (3.4)  and  r 

n 


(4.17) 


(H(J.t0.Vh)«)i  -  — 


/  (H+(ui,tN,h)(a)I  -  J+)_1*  ) iduj 


/  <H_(-X,t0,h)UI  -  j')_1z-)idl 


-19- 


where  r+  <=  {z  8  C|Rez  >0},  r_  c  {z  8  C|Rez  <  0}  holds.  Since 


(4.18) 


max  Mol  -  J+)  I  ,  maxl(XI  -  J  )  '  I  <  const. 


uer. 


xer 


holds,  the  estimates  given  in  Chapter  3  can  be  used  because  they  were  formulated  uniformly 
for  -X,(D  in  compact  subsets  of  the  left  half  plane. 

By  evaluating  (4.12)  at  the  boundaries  (i  »  0  resp  i  -  N)  we  get  the  block  system 


r 

0 

BE 

BE 

o 

j" 

_e 

(4.19) 

r~ 

S(T)E[  *  ]  S(T)E 


*0,K-1(-J  *h> 


.-1. 


6  -  BE(H(J,tQ,tN,h)E  f)Q 


.-1. 


y(T)  -  S(T)E(H(J,t0,tH,h)E  f)N 


We  assume  that  (2.1),  (2.2),  (2.3)  with  A(T)  i  A  has  a  unique  solution  for  all  f'3  e 
r 


C(  [1,-J ),  8  8  R  •  Therefore  (2.15)  has  to  hold  which  implies  that  be 


is  non¬ 


singular.  From  (2.30)  we  conclude  that  S(T)e[  *  ]  has  a  bounded  inverse.  From  (4.16)  we 
conclude  that 

(4.20) 
and 

(4.21) 


IY  +  (J+,h)l  <  const,  max  |Y  ,(w,h)|  <  const. 

O.N-1  aej.  °'N_1 


,Yn  u  i<_J  <h>'  <  const,  max  |Y_  , (-X,h)|  . 

0.N-1  xer  0,N-1 


Now  let  X  be  the  eigenvalue  of  j”  which  is  nearest  to  the  imaginary  axis,  such  that 


ReX  ■  -XmXn  and  take  r_  such  that  for  some  small  8  >  0 


(4.22) 


dist(r  ,X)  ”  6  and  dist( r  ,  {z  8  C|Rez  «  0})  »  X  -  5 

—  min 


.N-1 


holds,  we  now  choose  h  -  (h , > .  .  such  that  h,  <  h(X  .  -S.c.t,)  defined  in  (3.21). 

i  i*u  1  min  1 
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Then 


with  A(T )  i  A  we  get 


The  number  of  necessary  steps  N(e)  when  equality  holds  in  (4.38),  (4.39)  is  given  by 
(3.23),  (3.24) 

(4.41)  N(e)  *  c  —  ,  e  ♦  0 

/i 

which  is  comparable  with  the  number  of  steps  a  constant  stepsize  algorithm  would  need  for 
the  Integration  of  a  constant  interval  problem  in  order  to  achieve  an  0(e)  accuracy. 

Of  course  the  second  term  in  the  error  estimate  (4.36)  can  be  reduced  by  adding 

1 

gridpoints  t  >  (— - c^~ — -  in  — )a+1  and  by  forming  stepsizes  h  according  to  (3.21)(b) 

1  min-®  c_  l 

(second  and  third  branch  of  h  ).  An  estimate  for  the  number  of  gridpoints  in  the  case 

hi  “  Vin-*'6'*!*  iB  <jiv®n  (3.30). 
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(5.11) 


(5u)  i*  V2 "  V2  G(tu  v2  1  <ut+i+ui) '  Gu  “  ( <5u)  *♦  V2  >  w 


has  been  set.  £+  6  C  and  e  C  hold. 

Prom  (4.17)  and  the  Lemmas  3.2  and  3.3  we  get 


max  l(H<J,t  ,t  ,h)Gu)  I  < 
£-I(1)H  1  N  * 


(5.12) 


He 


<  const. 161 


It  t  1  <1  +  2h£tl+1A)  ,u£* 

lW  £-I(1)(N-.1)  Z  *  *  /2  t-KDN 


_1_ 
1  ^  ottl 


choose  T  «  t^  ■  (— - in  — )  and  h^  <  ''(X^^-S*  defined  in  (3.21). 

min 


X  >0  is  again  the  modulus  of  the  real  part  of  that  eigenvalue  of  A(  <■>)  which  is 
min 

closest  to  the  Imaginary  axis  of  all  eigenvalues  of  »( “)  with  negative  real  part  and 
6  >  0  is  sufficiently  small. 

Since  h(  X1  ,e,t4+1 )  >  hd^e.tj)  we  derive  from  (3.21)  that  h(  X^^  -  S.e.t^) 

<  T_a  such  that  1  +  y  httj+1/  *  const  holds.  He  define  the  operator 
(5.13)  HJh.tj.t^)  t  C  ♦  C 


......  n(N-I )  ,  .  N- 1 

such  that  for  x  e  ,  x  -  (xt+^  )t-J 


(5.14) 


HCh.tj.tjjJx  - 


(H(J,tI.tN.h)x)I 


OKJ.tj.t^.hJx^ 


[vv 


holds.  From  '5.12)  we  get 

(5.15)  IH(h, t  i t  )G 1  <  const.  161. 

I  N  l 

where  1*1  denotes  the  max-norm  for  vectors  in  the  respective  resp.  the  associated 

matrix  norm.  Therefore  the  operator 

(5.16) 

is  invertible  for  tj  <  tN  sufficiently  large.  He  define 


I  -  H(h,t.t  )G  :  c",N-I+1>  ♦c”'"-1*1’ 
I  N 
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(5.17)  (a)  Y*  U(J%)  - 

X  #  N 


such  that 


C-i<j+'h) 


Cb)  Y  „<-J  ,h)  - 

X#N 


'h,e 


'h,e 


1+1 

• 

Vi 

% 


~  -i  +  ,_+ 


~  -i  - 


-  (I-H(h,tj,t^)G)  »h)5+  +  (I-fKh.tj.t^G)  Yi#n(-J  .h)t_  + 


+  (I-H(h,t1,tN)G)"1HE'1f 


holds.  In  order  to  obtain  ^  the  difference  equation  (5.7)  has  to  be  Bolved 

backwards  with  qlven  . 

^1  ot  “1 

Therefore  h0',**'hI_1  have  to  be  chosen  such  that  (I  -  t“+  A(tt+  y) )  exist 

for  t  »  0 ( 1 ) ( I—  1 ) •  From  this  and  (5.18)  we  get 

_+ 


(5.19) 


U1  "  Zi  5+  +  Zie-  +  *i(f>  ' 


1  -  0 (  1  )N 


where  Z*  is  a  n  x  r+  ,  Z^  is  a  n  x  r_  matrix  and  r^lf)  0  c"  . 
The  block  system 


(5.20) 


BEZ0 

BEZ0 

f'+)- 

B  -  BEZQ(f) 

S(T)EZ* 

S(T)EZ^ 

1  5.  J 

y(T)  -  S ( T ) EZ  ( f ) 

N 

results  by  evaluating  at  the  boundaries  1  «  0  resp  1  •  N  and  by  using  (5.8),  (5.9).  We 
set 


(5.21)  (a)  Z 

From  (5.18)  we  get 

(5.22) 

and  since  Ht(h,t  ,^>1 


I/N 


z+ 

z“ 

I 

I 

• 

_ 

• 

• 

,  (b)  ZT  u  - 

• 

• 

I/N 

• 

+ 

ZvT 

N 

N 

ZI  ,N*  l  (H(h'VVG~>  ^,N(J+'h) 


t-0 

«  const  as  e  ♦  0 
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,ZI,W  -  <«(J+'h>l 


(5.23) 


<  const.  max  IG(t  y)\ 


t-KIHN-l) 


p 

—  « 

r+  + 

r+  + 

1 

Vi,N-1,J  'h) 

♦ 

Vi,s-i(J  'h) 

0 

0 

-c(T-t ) 

<  const.  max  (lG(t  y  ) la  )  . 


This  follows  from  (4.12)  and  Lemma  3.4.  It  is  easy  to  check  that  the  right  hand  aide  of 

(5.23)  converges  to  zero  as  T  ♦  •  (s  *  0).  Therefore 

(5.24)  ZN  “  f  0  1  +  0(1)  '  6+0 

follows  and 

S(T)SZ*  -  S(T)B[  q  ]  +  0(1) 

is  nonsingular  for  e  sufficiently  small  because  of  (2.30)  and  the  matrix  in  the  (1.1) 
position  of  (5.20)  is  bounded  as  e  ♦  0  . 
we  conclude  from  (5.18) 


(2.25) 

In  Chapter  2, 


Zi,N  -  H<h'VV«VH  +  • 

(2.16)  it  was  noted  that  the  general  solution  of  the  homogenous  problem 

u’»  t°(J+G(t))u 
u  8  CC  [1,H  ) 


can  be  written  as 

(5.26)  *_<t>  -  ((I  -  HgTV  •)[  z°  ])(t)  ,  t  >  y 

r 


where  the  solution  operator  H  is  defined  in  (2e11)  and  $(t)  is  as  in  (2e9).  We  get 

(5.27)  tjijt)  -  (HGi(p_)  (t)  +  *(t)[  °  ]  ,  t  >  i 

r 

and  define  the  vector  ^  for  y  m  t 

I  /  N  I 


(5.28) 


r  o 

♦_(tx) 

(HGi)i_){tI) 

♦‘V  I 

r 

• 

e 

• 

e 

at 

• 

4- 

• 

• 

e 

• 

*-(V 

(HG*_)(tN) 

•>v  k 

J  V 

Lr-J 

(HG*->I,N  ♦l.N 
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By  Proceeding  es  in  Chapter  4  we  get 


Theorem  5 ■  1 .  Assume  that  A  8  C((1,*])i  A*,  A"  8  C([1,»))  n  IMM,"])  and  that  f 
fulfills  (2.13)  with  P,  P' ,  P"  e  C(  (1/*))  n  L^( [  1 , «))  •  Let  for  some  e  sufficiently 


snail 
(5.36  ) 

hold  with  some  fixed  small  6 
(5.37)  h£ 


T(e)  -  t, 


(T 


a+1 


X  .  -«) 

min 


in  -) 
E  1 


1_ 
1  \0tM 


and  assume  that 


<  r0/e 


for 


fci  4  * 


c0  >  0 


(5.38)  h±  <  /e  t~aexp(-  -  t^*"1 )  ,  y  <  tt  <  T(e) 

hold  for  some  fixed  y  sufficiently  large.  Then  if  the  matrix  (2.18)  is  nonsingular  the 
Box-scheme  (5.1),  (5.2),  (5.3)  is  uniquely  soluble  for  e  sufficiently  small  and 

(5.39)  max  lx  -  y(t  ) I  -  0(e) 

i-0(1)N 


holds  for  y(T)  r  0  . 

If  equality  holds  in  (5.37), 
(5.40) 


(5.38)  the  number  of  steps 
n(e)  -c,^-,  e*0  . 


N  -  N(E)  fulfills 


The  condition  number  of  a  nonsingular  A  is  defined  by 

(5.41)  x(A)  -  IAIIA-1!  . 


Then  the  condition  number  of  the  difference  operator  1^  (given  by  (5.4),  (5.5),  (5.6) 
fulfills  the  estimate 


(5.42)  x(L.  )  <  const  — —  ~  const.  N(e) 

h  A 

if  equality  holds  in  (5.37),  (5.38). 

This  holds  since  II.  I  <  const.  — —  and  because  of  the  stability  estimate  (5.35). 

^  ft 

(5.42)  is  a  very  moderate  condition  number  and  therefore  (5.4),  (5.5),  (5.6)  can  be 
safely  solved  (by  partial  pivoting  using  S0LVEBL0CK(de  Boor  and  Weiss  (1980)). 


a  >  0 


6.  Nonlinear  Problems. 

Me  consider  the  'infinite'  problem 

(6.1)  y'  -  taf(t,y).  1  <  t  <  -  , 

(6.2)  b(y( 1 ) l  =  0 

(6.3)  y  S  C(  [1,«]  ) 
and  the  Box-scheme 

X  “X . 

(6#4,  -  *r+  v2  f  (tit  v2'  v2  (xi+i+xi)  > 

(6.5)  b(xQ)  “  0 

(6.6)  S(T)x„  =  S(T)y*(“) 

where  T  «  tj,  holds.  The  asymptotic  boundary  condition  S(T)  is  considered  with  regard 
to  Chapter  2. 

As  mentioned  in  Chapter  2 

(6.7)  f(“,y(<»)}  «  0 


0(1)(N-1) 


has  to  hold.  He  now  assume  that  there  is  an  isolated  zero  y  (»)  and  that  f  (°°,y  (<*>)) 


has  the  Jordan  form  J  obtained  by 

(6.8) 

where  J  fulfills  (2.5).  Moreover  we  assume 
(6.9) 


fy(«o,y  (oo)  )  m  EJE 


-1 


1 


)»  U  >  X 


min 


2  *  *  rrM 

f  e  C^(Cx(1,y  (-)),  f(t,y  (-))  -  0{e 
r_ 

(6.10)  b:  R  R  ;  b,by  are  locally  Lipschitz  continuous  in  Rn 


and  that  the  problem  (6.1),  (6.2),  (6.3)  has  an  isolated  solution  y  (t)  ♦  y  (~)  as 
t  ~  .  The  isolatedness  means  that  the  linearized  problem 


fy(t,y  (t))Z 


by(y  (l))z(l>  -  0 


z  e  c<  [1,-1 ) 

has  only  the  trivial  solution  z  =  0.  Then  we  conclude  from  de  Hoog  and  Weiss  (1980a)  that 
the  approximating  problem 

(6.11)  x'  -  t“f(t,x) 

(6.12)  b(x( 1 ) )  -  0 


(6.13) 


S(T)x(T)  -  S(T)y  (®) 


with  (2.29) ,  (2.30)  are  locally  (around  y  ( t ) )  uniquely  soluble  for  T  sufficiently  large 


and 

(6.14)  lx  -  ylr,  .  <  const.  IS ( T) (y(T)  -  y  (»))l 

L  I 

* 

holes.  An  estimate  for  S(T)(y(T)  -y  («))  can  be  obtained  from  (2.21).  Possible  choices 
for  S(T)  are  discussed  in  Lentini  and  Keller  (1980). 


1 

Now  we  choose  T  -  (-■■■—- — r—  in  — 1  for  e  sufficiently  small  and  apply  the 

min  6 

nonlinear  stability  theory  given  in  Keller  (1975)  with  e  as  a  parameter.  The  result  then 
follows  from  the  stability  estimate  (5.35)  for  linear  problems  and  we  merely  state  it  in 
Theorem  6. 1 .  Under  the  given  assumption  the  Box-scheme  is  convergent  for  e  ♦  0  to  the 
locally  unique  solution  of  (6.4),  (6.5),  (6.6)  if  stepsise  sequence  h^  fulfilling 

(6.15)  h^  <  cfl/e  ,  tt  <  y 


(6.16)  ht  <  /l  t~aexp(  Y  <  <  T (6) 

is  chosen.  The  estimate 


(6.17)  max  ly(t  )  -  x  I  •  Ole)  ,  e  ♦  0 

i-0(1)N 

holds.  The  Newton  procedure  for  (6.4),  (6.5),  (6.6)  is  quadratically  convergent  for 
starting  values  in  a  sphere 


K 


{x8Cn(N+1,|.x- 


y<V 

y<V 


i  < 


where  £  is  independent  of  e  . 
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ff'TT* - IV 


■  ,'Ktam m 


■  ifiaiiS ...... 


7.  Algorithm  1 

A  difficulty  that  might  arise  is  that  T(e)  -  (  —  - - —  In  —  is  rather  small 

-4  -8  "ln  '  -8 

for  reasonable  e's  (10  <  e  <  10  ).  For  example  assume  a  -  1,  "  2»  e  «  10 

Then  T  «  4. 

If  y(t)  has  not  'reached'  its  asymptotic  state  at  T  <•  4  no  good  approximation  by 
the  finite  interval  problem  posed  on  [1,4]  can  be  expected.  However  a  significant  increase 
of  T  (such  that  y(T)  is  reasonably  close  to  y(«) )  which  corresponds  to  an  enourmous 
decrease  of  e  would  imply  a  large  increase  in  the  amount  of  labor  since  N(e)  * 

-  0(  e”  1/2  ). 

A  reasonable  way  to  overcome  this  difficulty  is  to  set 


1 

(7.1)  T(e)  -  (tt — —  in  -  )at1 

(Xmin  ~  5)  e 

where 

X  -  « 

(7.2)  K  -  max  |(y(t)  -  y  ( -) )  exp  -  t0^  J| 

te[1,«]  ^ 

has  been  set. 

Then 

(7.3)  ly(T(c))  -  y(«)l  <  e  for  all  e  >  0  . 

The  meshsize  sequence  (6.15),  (6.16)  can  still  be  used  on  the  whole  interval  [ 1 ,T( e) ] . 
The  error  estimate  (6.17)  takes  the  form 

(7.4)  max  lx  -  y(t  )l  »  0(ke),  e  ♦  0 

1-0(1  )N 

-Vo 

and  N( e)  would  still  be  an  0(e  2  ). 

In  the  case  that  the  function  f  ,  which  sets  up  the  differential  equation,  is 
independent  of  t  ,  the  choice 

(7.5)  S(T )  =  S  -  [Ir  ,0]E_1 
(and  y(T)  =  y  -  Sy(»))  implies 

(7.6)  lx  -  ylr,  <  const  ly(T)  -  v(»)l2 


(see  Lentini  and  Keller  (1980)).  Therefore  we  can  choose 


(7.7) 


t 

TUl.^t^rrln^r 

min 

and  (7.4)  still  holds. 

X  way  to  construct  an  adaptive  code  baaed  on  the  given  theory  is  the  following. 

At  first  choose  a  T  so  large  that  y(T)  is  reasonably  close  to  y(  <•) .  This  might  be 

done  by  using  physical  information  on  the  solution  of  the  infinite  problem.  Then  choose  a 

coarse  grid  on  (1.T)  such  that  the  meahsizes  increase  exponentially  as  t  +  T  .  The 

solution  of  the  equations  (6.4).  (6.5).  (6.6)  gives  an  initial  guess  (x  .••♦.x  ).  Now  we 

0  N 

calculate 

\  -  6 

(7.8)  ic  ■  max  |(x  -  y(  «■)  )exp(— -  t,)|  . 

1-0(1  )N  on  l 

Here  5  should  be  an  input  parameter.  Then  we  set  y  -  tj  (for  the  meaning  of  y  see 
Theorem  6. 1 )  such  that 

(7.9)  max  |x  -  y(«)|  <  ~  t  |x  -  y(»)|  >  —  for  some  j  <  I 

*-I(1)N  *  m 

where  H  is  also  an  input  parameter. 

For  i  »  1(1) (1-1)  we  calculate  the  first  order  term  of  the  local  error  JL  using 
the  x^'s  as  bentlni  and  Pereyra  (1977)  did.  On  [1,tj]  we  define  an  equidlstributing 
mesh  proceeding  as  in  the  just  cited  reference  given  an  e  >  0  .  On  [tj,  T(e)l  where 
T( e)  is  given  by  (7.1)  resp.  (7.7)  with  k  set  for  k  .  we  choose  the  meshsizes  h^  as 

h  -  /e  t  °  exp(  •  ”  — -  t?^1)  ,  i  -  l( 1 ) (N-2) 

1  l  v  or*- 1  i 

Vi  “  ?(e)  •  Vi 

and  solve  (6.4),  (6.5),  (6.6)  on  (1,T(e)l  using  the  constructed  mesh. 

A  standard  error  estimation  algorithm  then  checks  whether  a  given  accuracy  has  been 
achieved.  If  yea  then  the  algorithm  stops,  if  not  we  calculate  from  the  just 

obtained  solutions  of  the  Box-scheme  and  calculate  a  new  grid  with  "  e/2  separately 

°"  [1,YNEW!  *  fYNEW  '  ?<ENBWn  ^  in  th®  flrSt  iteration- 

This  interative  procedure  stops  as  soon  as  the  required  accuracy  is  obtained. 

Numerical  experiments  will  be  reported  in  a  subsequent  paper  solely  concerned  with 

computational  aspects  of  'infinite'  boundary  value  problems. 
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